Linear perturbative theory of the discrete cosmological N-body problem 
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Abstract 

We present a perturbative treatment of the evolution under their mutual self-gravity of particles 
displaced off an infinite perfect lattice, both for a static space and for a homogeneously expanding 
space as in cosmological N-body simulations. The treatment, analogous to that of perturbations to a 
crystal in solid state physics, can be seen as a discrete (i.e. particle) generalization of the perturbative 
solution in the Lagrangian formalism of a self-gravitating fluid. Working to linear order, we show 
explicitly that this fluid evolution is recovered in the limit that the initial perturbations are restricted 
to modes of wavelength much larger than the lattice spacing. The full spectrum of eigenvalues of 
the simple cubic lattice contains both oscillatory modes and unstable modes which grow slightly 
faster than in the fluid limit. A detailed comparison of our perturbative treatment, at linear order, 
with full numerical simulations is presented, for two very different classes of initial perturbation 
spectra. We find that the range of validity is similar to that of the perturbative fluid approximation 
(i.e. up to close to "shell-crossing"), but that the accuracy in tracing the evolution is superior. The 
formalism provides a powerful tool to systematically calculate discreteness effects at early times in 
cosmological N-body simulations. 

PACS numbers: 98.80.-k, 05.70.-a, 02.50.-r, 05.40.-a 



I. INTRODUCTION 



The standard paradigm for formation of large scale 
structure in the universe is based on the growth of 
small initial density fluctuations in a homogeneous and 
isotropic medium (see e.g. 0). In the currently most 
popular cosmological models, a dominant fraction (more 
than 80 %) of the clustering matter in the universe is as- 
sumed to be in the form of microscopic particles which in- 
teract essentially only by their self-gravity. At the macro- 
scopic scales of interest in cosmology the evolution of the 
distribution of this matter is then very well described by 
the Vlasov or "coUisionless Boltzmann" equations cou- 
pled with the Poisson equation (see e.g. 0). A full 
solution, either analytical or numerical, of these equa- 
tions starting from appropriate initial conditions (IC) is 
not feasible. There are, on the one hand, various pertur- 
bative approaches to their solution (for reviews see e.g. 
H H), which allow one to understand the evolution in 
some limited range (essentially of small to moderate am- 
plitude fluctuations). On the other hand, there are cos- 
mological N-body simulations (for reviews see 



which solve numerically for the evolution of a system of 
N particles interacting purely through gravity, with a 
softening at very small scales. The number of particles 
N in the very largest current simulations Q is 10^", 
many more than two decades ago, but still many orders 
of magnitude fewer than the number of real dark matter 
particles (~ 10^° in a comparable volume for a typical 
candidate). The question inevitably arises of the accu- 
racy with which these "macro-particles" trace the desired 
correlation properties of the theoretical models. This is 
the problem of discreteness in cosmological N-body simu- 
lations. It is an issue which is of considerable importance 
as cosmology requires ever more precise predictions for 
its models for comparison with observations. 

Up to now the primary approach to the study of dis- 
creteness in N-body simulations has been through nu- 
merical studies of convergence (see e.g. 0, i-^-j 
one changes the number of particles in a simulation and 
studies the stability of the measured quantities. Where 
results seem fairly stable, they are assumed to have con- 
verged to the continuum limit. While this is a coherent 
approach, it is far from conclusive as, beyond the range 
of perturbation theory, we have no theoretical "bench- 
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marks" to compare with. Nor is there any systematic 
theory of discreteness effects, e.g., we have no theoreti- 
cal knowledge of the A'' dependence of the convergence. 
Given that typically N is varied over a very modest range 
(typically one or two orders of magnitude) compared to 
that separating the simulation from the model (typically 
70 orders of magnitude) there is much room for error. 

Different mechanisms by which discreteness effects 
may make the evolution of N-body simulations different 
to that of the fluid limit have been discussed in the lit- 
erature. A very basic consideration is that of the dis- 
creteness effects introduced already in the IC, before any 
dynamical evolution. Indeed there are necessarily dis- 
crepancies between the correlation properties of the dis- 
cretised IC and those of the input theoretical model, as 
there are intrinsic fluctuations associated with the par- 
ticles themselves. For analysis and discussion of these 
effects see, e.g., ful H ES Q E 113 ■ What is probably 
the most obvious effect of discreteness, and certainly the 
one most emphasised in the literature, is two body coUi- 
sionality: pairs of particles can have strong interactions 
with one another, which is an effect absent in the coUi- 
sionless limit. For a naly sis and discussion of these effects 
see, e.g., [H El IH EE El ■ 

In this paper we present an approach which allows in 
principle a systematic understanding of discreteness ef- 
fects between these two regimes, in the evolution from the 
IC up to the time when two body collisions start to oc- 
cur. We do so by developing a perturbative solution to 
the fully discrete cosmological N-body problem, which 
is valid in this regime. This essentially analytic solu- 
tion can be compared to the analogous fluid {N — s- oo) 
solution, and one can understand exhaustively the modi- 
fications introduced, at a given time and length scale, by 
the finiteness of N. While the usefulness of the approach 
is restricted to the regime of validity of this perturbative 
approach, we can gain considerable insights into the ef- 
fects of discreteness and how they introduce error. Some 
of the essential results have already been briefly reported 
in In this paper we describe in much greater detail 
the perturbative method used to describe the evolution, 
and evaluate its regime of validity by extensive compari- 
son with numerical simulations. In a forthcoming paper 
[23j | we will discuss the application of this method to 
the study of discreteness effects in N-body simulations, 
providing precise quantifications of these effects in the 
regime in which our treatment is valid. 

The perturbative scheme we employ is one which is 
well known and standard in solid state physics, as it is 
that used in the treatment of perturbations of a crystal 
about the local or global minimum of its internal energy. 
Indeed the class of cosmological N-body simulations we 
consider are those which start by making very small per- 
turbations to particles initially placed on a perfect simple 
lattice ISISIllEllllEl. Up to an overall change in 
sign, our perturbative scheme is precisely that one would 
use for the analysis of the extensivel y st udied Coulomb 
lattice or Wigner crystal (see e.g. H^,!!!!), of iV particles 



on a lattice interacting by a pure unscreened Coulomb 
force. At linear order (harmonic analysis) one simply 
solves a 3A^ x 3A^ eigenvalue problem to determine the 
eigenmodes and eigenvalues of the displacements off the 
crystal. This can be done at very low cost in computa- 
tional resources because of the symmetries of the lattice. 
Stable (i.e. dynamically oscillating) modes in one prob- 
lem become unstable (growing and decaying) modes in 
the other problem, and vice versa. One consequence of 
this which we will discuss briefly here, and more exten- 
sively in |23|, is that, for what concerns discreteness in 
N-body simulations, there are qualitatively different fea- 
tures on the simple cubic (sc) lattice compared to the 
body centered cubic (bcc) or face centered cubic (fee) 
lattice. 

A crucial step in our analysis is evidently the compari- 
son of our solutions for the evolution with those obtained 
from the treatment of the continuous self-gravitating sys- 
tem. This latter problem can also be solved in a limited 
range with a perturbative treatment of the equations of 
a self-gravitating fluid, which are obtained by trunca- 
tion of the coUisionless Boltzmann equation. Given that 
our perturbation scheme works with the displacements of 
the particles, one might anticipate that the appropriate 
perturbation scheme to compare with is that given in the 
Lagrangian formulation of these fluid equations, in which 
the evolution is described in terms of the trajectories of 
fluid elements [S^. We show explicitly, at linear order, 
that this is the case: taking the limit in which the per- 
turbations to the lattice are of wavelengths much larger 
than the lattice spacing £ the evolution described by our 
scheme maps precisely on to that at the same order in the 
Lagrangian description of the fluid. The Zeldovich ap- 
proximation, which is simply the asymptotic form of this 
solution, can then be understood in very simple analogy 
with the long-wavelength coherent "plasma oscillations" 
in a unscreened charged plasma. 

The paper is organized as follows. In the next section 
we describe the perturbative treatment for perturbations 
off a perfect lattice, for the specific case of gravity. We 
work firstly, for simplicity, in a static Euclidean universe, 
giving the explicit expressions for the evolution from gen- 
eral IC (i.e. any perturbation from the lattice). In the 
following section we explicitly solve the evolution for the 
case of a simple cubic lattice, and discuss in detail the 
structure of the spectrum of eigenvalues. In Sect. IV we 
generalize our results to the case of an expanding uni- 
verse, and then show explicitly the recovery of the fluid 
limit given by the solution at the same order of the La- 
grangian formulation of the equations of a self-gravitating 
fluid. In the next section we present a comparison of the 
evolution described by our approximation with that of 
full numerical N-body simulations. We consider both un- 
correlated initial perturbations (a "shuffled lattice" ) and 
a set of highly correlated perturbations (with a power 
spectrum of density fluctuations fc~^), like that in cur- 
rent cosmological models. We verify that the agreement 
is very good, and that the evolution is traced with consid- 
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erably better accuracy than by the fluid hmit at the same 
order. Both approximations break down when particles 
start to approach one another (i.e. "shell-crossing" in 
fluid language), which we parametrise through an appro- 
priate statistical measure. In the final section we sum- 
marise our results and discuss various further develop- 
ments of the method presented here which could be pur- 
sued, notably the extension of the perturbation scheme 
to higher than linear order. This may throw further light 
on how insights about discreteness gained using this for- 
malism, which will be discussed at length in l2i|| , extend 
into the regime of highly non-linear evolution. We also 
discuss briefly the possible interest of solving the cos- 
mological N-body problem on the bcc lattice, as well as 
the possible extension of our method to IC generated on 
"glassy" conflgurations. 



II. LINEARIZATION OF GRAVITY ON A 
PERTURBED LATTICE 

In this section we start by defining and studying some 
general properties of the gravitational potential and force 
of an infinite system of point particles. We then consider 
the particular case of a perturbed infinite lattice in a 
static Euclidean space, the generalization to an expand- 
ing universe being given in Sect. IIVI Since the force is 
zero in the unperturbed lattice, the dominant contribu- 
tion to the force in the perturbed case is linear in the 
relative displacements of the particles. In the last sub- 
section, we consider the equations of motion resulting 
from this linearized force. 



A. Definition of tlie force and tiie potential 

Let us consider carefully first the definition of the grav- 
itational force in an infinite system of point particles of 
equal mass m. We will assume that this system (either 
stochastic or deterministic) is characterized by a well de- 
fined mean number density no > 0, and mass density 
Po — rnuQ. The gravitational potential of a particle, per 
unit mass, at r, due to the particles in a finite volume V, 
is: 



(r) = -Gm ^ ■ 



1 



(1) 



where the sum is over all the particles contained in the 
system, and V{V, r) is the window function for the volume 
V, i.e.. 



V{V,r) 



1 ifreV, 
0, otherwise. 



(2) 



The force per unit of mass (i.e. the acceleration), due 
to these same particles, is given by the gradient of the 
potential: 



Taking the infinite volume limit V oo, neither the 
gravitational potential nor the gravitational force 
©, are well defined. In the first case the result diverges, 
while in the second it may be finite or infinite, but its 
value depends on how the limit is taken ^ . 

In Euclidean spacetinie this behaviour in the infinite 
volume limit may be regulated by the introduction of a 
negative background — the so-called Jeans swindle (see 
e.g. 0,113) — that the potential is defined as 



(r) -G lim 



-Po 



(4) 



This modifies the usual Poisson equation to 

V20(r)=4^G(p(r)-po). (5) 

The expression Q is well defined^, provided (i) that 
the limit y ^ oo is taken in a physically reasonable 
way'^, and (ii) that the fiuctuations in the system are 
sufficiently rapidly decaying at large scales . In the cos- 
mological context this negative background appears nat- 
urally as a consequence of the expansion of the universe 
(see Sect.lTvTl. 

The simulations of self-gravitating systems we are in- 
terested in are performed using a finite cubic simulation 
box of side L and volume Vb = L^, subject to peri- 
odic boundary conditions. The force on a particle is thus 
computed not only from all the other particles inside the 
simulation box, but also from all the copies of the par- 
ticles contained in the "replicas". The reason for using 
these boundary conditions is that they introduce the in- 
evitable finite size effects without breaking translational 
invariance: every particle can be considered to be at the 
centre of the finite box and therefore sees the boundary 
in the same way. The infinite system we consider is thus 
an infinite number of replicas of a finite cubic box, with 
a negative background as described above to make the 
force well defined ^ . In this case the gravitational poten- 
tial may be written as 



(r) = lim [0h(r) -I- cj)p{v)\ , 



(6) 



r(r) = -V0(r) 



(3) 



^ F(r) is a conditionally convergent series. 

^ For a more detailed discussion of the gravitational force in infinite 
systems see also j31j. 

^ E.g., taking the infinite volume limit in compact sets. 

* If P(k) is the power spectrum of density fluctuations, it is sim- 
ple to show, using the modified Poisson equation Eq. that 
convergence of the fiuctuations in the gravitational potential re- 
quires lim(;_,o = for n > 1. For finite fiuctuations in 
the force one requires n > — 1. 

^ Note also that, because the system is just a lattice when con- 
sidered at scales larger than the box size, the fluctuations are 
always sufficiently suppressed at large scales so that the gravita- 
tional force is well defined. Thus any possible divergence in the 
fluctuations of force will be regulated by the box size L. 
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where 

0fc(r) = Gpo / dV^-V(l/,r') 
is the contribution from the background, and 



Note that the contribution from the background (|10(l is 
identically zero if one takes a window function with in- 
(7) version symmetry in r (e.g. a sphere or cube centred on 
r). 



(8) B. Linearization of the gravitational force 



the contribution from the particles. Here the sum over r' 
is restricted to the particles in the box, while the other 
sum, over the three integers n (i.e. over the images of 
r'), has a "*" to indicate that the term r' = r is excluded 
when n = 0. 

The gravitational force is: 



F(r) = ^lim [Ff,(r)+Fp(r)], 



(9) 



We consider the infinite lattice generated by the repli- 
cation of a sc lattice of volume Vb of side L with N ele- 
ments, i.e., whose lattice vectors are R = (mi,m2,TO3)£ 
with m, e [0, N^^^ - 1] n N and ^ = L/N^^^ is the lattice 
spacing^. This lattice (with a particle at each site) is 
now perturbed by applying displacements u(R) to each 
particle R, so that the new positions of the particles can 
be written as 



where 



and 



F.(r).Gpoy^/V^-^V(F,r') 



Fp(r) 



-Gm 



nL 



(10) 



V{V,r' + nL). (11) 



The generalization^c ons pres ented here to any 

Bravais lattice is straightforward (see e.g. l32l t . 



r = R + u(R) 



(12) 



The "particle" term in the gravitational force [i.e. 
Eq. (|ll|l ] can then be expanded order by order in Taylor 
series about its value in the unperturbed lattice. At lin- 
ear order in the relative displacements u(R) — u(R') we 
obtain 



Fp(r) = -Gm ^ I 



R - R' + nL u(R) - u(R') [u(R) - u(R')] • [R - R' + nL] 



n,R 

X V(V",R' + nL 



R-R' + nL|3 |R-R' + nL|3 



|R- R' + nL|' 



(R - R' + nL) 



(13) 



The first term in this sum 

n,R' ' 

has the poor infinite volume behaviour which is regu- 
lated, as discussed above, by the contribution coming 
from the background Eq. H1U|I . The total linearized force 
is then also well defined, and given by the infinite vol- 
ume limit of Eq. H13|l summed with Eq. (|10|) . In the 
case that we choose to calculate using the infinite vol- 
ume limit of a volume V with inversion symmetry in r 
(i.e. the displaced position of the particle), the full lin- 
earized force is thus given by Eq. (|13|l . If, however, we 
choose to sum in a volume with inversion symmetry in 
the lattice site R, it is simple to show that Eq. (|14|l is 
identically zero. The background term then contributes, 
with the sum [ H1U|) + (|14|l ] remaining invariant. 



The convergence criterion for each term of (|13|) is 

|R-R'| > |u(R) -u(R')|. (15) 

Note that the validity of the power expansion does not 
depend on the displacement of the particle R on which 
we compute the force, but on relative displacements of 
the particles at the position R and R'. Under the ac- 
tion of the gravitational interaction, the displacements 
u(R) will typically grow so that the condition Eq. H15|) 
is violated after some time. However when some pairs 
of particles no longer satisfy condition (|15() . it may nev- 
ertheless continue to apply for the rest of the particles 
and (|13|) may remain a sufhciently good approximation 
to the force. In order to have a precise characterization 
of the regime of validity of the approximation applied to 
follow the dynamical evolution of a perturbed lattice, it 
is necessary to compare the results with those obtained 
from evolution under full gravity. We will perform such 



5 



a comparison in Sect. El using N-body simulations. 

It is convenient to write the linearized force just dis- 
cussed in terms of the so-called dynamical matrix 2?(R) 
(see e.g. [Uli^): 



F(r) = 5^I?(R-R')u(R') 



(16) 



R' 



This matrix has the following properties: it is a com- 
plete symmetric operator, i.e., I?^^(R) = 'D^^{—'R) with 
inversion symmetry, i.e., I?^i/(R) = I?p^(— R). Further, 
since the same displacement applied to all the particles 
produces no net force, we have X^r^m"!-^) = 0. For 
any pair interaction potential v(r)it is straighforward to 
show that it can be written as 132. 1331 




FIG. 1: Computation of the diagonal terms of the dynamical 
matrix at R = 0. 



V^,(R ^ 0) = <9^a,w(R) 
V^,{R - 0) = - ^ d^dM^') 



R'#0 



where 



w{r) 



drf^dr^ 



and u;(r) is the periodic function defined as 
w{v) = ^ i;(r + nL), 



(17a) 
(17b) 



(18) 



(19) 



i.e., the potential due to a single particle and all its copies. 
For gravity we have ?;(r) = —Gm/r and Eq. (|19|l is im- 
plicitly understood to be regulated as discussed at length 
above by the addition of a uniform negative background. 
We will describe below, and in App. ^ how we use the 
well-known Ewald summation technique to explicitly per- 
form this sum. 

Equation (|17b|) gives the force on a particle, at first 
order in the displacements, when it is displaced and all 
the others remain unperturbed (see Fig. For gravity 
it is straightforward [SJ to show that 



47r 



(20) 



i.e., the linearized force fs(r) on a particle due only to its 
own displacement u with respect to the rest of the lattice 
is 



f.(r) 



47r 



Gpou(R). 



(21) 



The simplest way to derive this result is by summing the 
force in spheres centred on the unperturbed position of 
the displaced particle. In this case it is straighforward 
to show, by symmetry, that the linearized direct particle 
contribution Eq. (|13|l is zero and the full force is given by 
the background term Eq. (|10|l . The result follows then 
simply from Gauss' law which gives that the force comes 
only from the region inside the sphere shown in Fig. ^ 



C. Equations of motion in a static Euclidean 
universe 

In this section we derive the equations of motion of 
the particles in the linear approximation, and then solve 
them. We treat first a static Euclidean space, giving the 
generalization to a cosmological expanding universe in 
Sect.lTvl 

Using Newton's second law and Eqs. ((T^ and ifTCI) we 
can write the equation of motion of the particles as: 



u(R, t) = XI " R')u(R', t), 



(22) 



R' 



where the double dot denotes a double derivative with 
respect to time. The expression H22II is a system of vec- 
torial coupled second order differential equations which 
can be reduced to an eigenvalue problem, using standard 
techniques. From Bloch's theorem |33| it follows that 
Eq. H22|l can be diagonalized by the following combina- 
tion of plane waves: 



1 



u(R,t) = -Xu(k,t)e'''^ (23) 

k 

where the sum over k is restricted to the first Brillouin 
zone, i.e., for a sc lattice to 



k^-„, 



(24) 



with n (ni,n2,n3) such that n, G [-N/2, N/2[r\Z. We 
denote by u(k, i) the Fourier transform of u(R, f): 



a(k,t) = Xu(R,t)e 



-jk R 



(25) 



R 



where the sum is restricted to the simulation box (i.e. 
without considering the replicas). Inserting Eq. H23|) in 
Eq. H22|) . we obtain for each k: 



{i(k,i) = P(k)u(k,t), 



(26) 



where 2?(k) is the FT of I?(R), defined analogously to 
(|25|l . From the properties of f (R) given above, it fol- 
lows that I'(k) is a real and symmetric operator which 
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satisfies 

47r 

\imV,,,{k)^—Gpo6^,. (27) 

k— »0 o 

We can now solve Eq. (|26() by diagonalizing the 3x3 
matrix 2?(k). For each k, this determines three orthonor- 
mal eigenvectors e„(k) with three associated eigenvalues 
(jj^j(k) {n = 1, 2, 3) satisfying the eigenvalue equation: 

P(k)e„(k) =c<.2(k)e„(k). (28) 

We can decompose each mode u(k, t) in the basis {e„(k)} 

as 

3 

u(k,t) = ^e„(k)/„(k,i). (29) 

ri=l 

Using Eqs. (j^ . and we get the following equa- 
tion for the coefficients /„(k, t): 

/„(k,t)=c^,2(k)/„(k,t). (30) 

Depending on the sign of oj^ (k) , we obtain two classes of 
solutions C/„(k, t) and T4(k, t). We choose them, without 
any loss of generality, satisfying 

C/„(k,to) = l, ;7„(k,to) = 0, (31a) 
K(k,to)-0, K(k,io) = l. (31b) 

The function J7„ (k, t) is associated with initial displace- 
ments and t4i(k, t) with initial velocities. If '^^^(k) > 
then 



where the matrix elements of the "evolution operators" 
V and Q are 

3 

7'^,(k,t) =^C/„(k,t)(e„(k))^(e„(k))„ (35a) 

n=l 
3 

Q^,(k,t) =^K(k,t)(e„(k))^(e„(k)),. (35b) 

n=l 

The operator V thus evolves the initial displacement field 
and Q the initial velocity field. 



III. DETERMINATION AND ANALYSIS OF 
THE SPECTRUM OF EIGENVALUES OF ©(k) 

In this section we describe the determination of the 
eigenvectors and spectrum of eigenvalues of the dynam- 
ical matrix for gravity. We then discuss the physical 
meaning of the results, notably identifying how the fluid 
limit is obtained and how corrections to this limit may be 
calculated. In this discussion we will use extensively the 
strict analogy between the case we are treating and the 
Coulomb lattice, or Wigner cystal, studied in condensed 
matter physics (see e.g. [l^). This is a system of posi- 
tively charged particles embedded in a negative neutral- 
izing background. The particles interact with a repulsive 
1/r potential instead of the attractive — 1/r potential of 
Newtonian gravity. Thus all our results are mapped onto 
those for the corresponding Coulomb lattice by making 
the formal substitution Gm^ — > — e^, where e is the elec- 
tronic charge*. 



UniKt) = cosh(w„(k)(t -to)), 

K (k, i) = sinh(tj„ (k) (t - to ) ) /^n (k) . 



(32a) 
(32b) 



If^2(k)<o 



C/„(k, t) = cos(Vl^2(k)|(t - to)), (33a) 
K(k,t) = sin(/RM(i - to))//RM- (33b) 

Whereas the modes (|32f) with positive eigenvalues cause 
an exponential growth of perturbation in the system, the 
modes (|33|l with negative eigenvalues leads to oscilla- 
tions. The evolution of the displacement field from any 
initial state u(R, to) is then given by the transformation 



u(R, t)^^Y. [^(^' *)"(^' *o) + 2(k, i)u(k, to) 



ik-R 



(34) 



But note that ©^^(k = 0) = Er^'mi^W = 0, i.e., V{k) is 
discontinuous at k = 0. 



A. Numerical computation of the spectrum of ©(k) 

The spectrum of the matrix 2?(k) must be computed 
numerically. The matrix P(Rj is constructed using the 
Ewald sum method [s^, IH, yj, to speed up the con- 
vergence of the sum. We continue to work here explicitly, 
as above, with a sc lattice of side L, with lattice spacing 
£ and N elements^. To determine the dynamical matrix 
we use the Ewald method to evaluate w{r) as given in 
Eq. H19|) , splitting it into two pieces using an appropriate 
damping function C: 



j(r) = ^ w(r + nL)C{\r + nL\,a) 

n 

^ w(r + nL)[l - C(|r + nL|, a)]. 



(36) 



° The potential we have used here for gravity has been defined per 
unit mass, i.e., in our notation v{r) = e'^/mr for the Coulomb 
lattice. 

® The generalization to a parallelepip ed b ox, and to other Bravais 
lattices, is straightforward (see e.g. |32| '). 



where a is a arbitrary "damping parameter" of which the 
result is independant. The function C(|r|, a) is chosen to 
be equal to unity at r = and rapidly decaying to zero 
as |r| goes to infinity. The first sum is then evaluated in 
real space and the second one in Fourier space, making 
use of the Parseval theorem (23, C being chosen so that 
the second term in Eq. is analytic at r = and thus 
rapidly convergent in Fourier space. A common choice 
for a l/r pair potential is 



C(|r|,a) = crfc(a|r|). 



(37) 



The expression for the function w is then: 

u;(r) = wW(r) (r). (38) 
In the gravitational case 

w^'^Ur) = -GmY ^erfc(a|r + nL|), (39a) 

y ' ^ \r + nL\ ^ ' " ^ ' 

n ' ' 

..(^)(r).-G^^g^exp(-g)cos[k.r], 

(39b) 



where Vb is the volume of the box and the wavevectors 
k are as in Eq. (|24|) . but with n ranging over all triple 
integers (i.e. not restricted to the first Brillouin zone). 
There is no k = term in the sum H39|l because of the 
presence of the negative background: when summed over 
all the particles, this term is equal to 



lim 00 (k) = - lim , 



(40) 



i.e., the k = mode of the potential (calculated from the 
Poisson equation in Fourier space) which is cancelled by 
the contribution coming from the negative background. 

The Ewald sum for the dynamical matrix can then be 
calculated directly using Eq. (flTjl and (p^ . The result, 
as in Eq. (|38|l . is divided in two parts: 



P(R) (R) +I?(''^'(R) 



(41) 



for which the explicit expressions are given in App. fXl 

For the results quoted here we have taken a — 2/L 
[33 ■ Using this numerical value of a, it is sufhcient to 
sum for 



|n| <3 



I. 1 Stt 
|k|<^. 



(42) 



to obtain a well converged determination of the dynam- 
ical matrix. The diagonalization calculation involves es- 
sentially N operations (where N is the number of parti- 
cles). It is perfectly feasible, with modest computer re- 
sources, to perform this diagonalisation for particle num- 
bers as large as those used in the largest current N-body 
simulations. 




FIG. 2: Spectrum of eigenvalues for simple cubic lattice with 
16^ particles. The lines correspond to chosen directions in k 
space. 



B. Analysis of the spectrum of eigenvalues in a 
simple cubic lattice 

We now describe the spectrum of eigenvalues of the 
dynamical matrix I?(R) for a sc lattice. As we have dis- 
cussed in the introduction, this is the lattice which is used 
very widely in N-body simulations of structure formation 
in cosmology. 

In Fig. 121 we plot the spectrum of a sc lattice, for N = 
16'^, obtained with the method outlined in the previous 
subsection. We show the normalized eigenvalues 



i(k) 



(43) 



as a function of the modulus of the k vectors, normalized 
to the Nyquist frequency fc^v — Tr/i- With this normal- 
isation the spectrum remains substantially the same as 
we increase the number of particles: the only change is 
that the eigenvalues become denser in the plot, filling out 
the approximate functional behaviours with more points. 
For our discussion here there is no interest in considering 
a greater number of points than that we have chosen. 

For each vector k there are three eigenvalues ij-'^(k), 
n — 1,2,3. Each family of eigenvalues (i.e. with same n) 
defines a surface, corresponding to the three branches of 
the frequency-wavevector dispersion relation. Sections of 
these surfaces are plotted for some chosen directions of 
the vector k in Fig. |21 



1. An expression for'D(k) and the Kohn sum rule 

Before proceeding further it is useful to derive some 
important results we will employ much in what follows. 
These are well known in the context of the application 
of this formalism in condensed matter physics (see e.g. 
[2^). First of all, we derive an analytical expression for 
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the dynamical matrix in Fourier space. Let us decompose 
in Fourier modes the function w{r) defined in Eq. (|19|l 

^W^l^E^W^'"'"' (44) 

where the sum over k is performed over all k space, i.e., 
not restricted to the first Brillouin zone and 

w{k) = [ dVu;(r)e-*''. (45) 
The derivatives of the periodic potential are 

w^Ar) = -^y] fc^fc.^'(k)e"^-"-. (46) 

k 

Using the definition of the dynamical matrix 

P^,(k) =^P^,(R)e-'^''-^ (47) 

R 

and Eqs. p7|l and (|46|l we obtain: 

^ k',R 

(48) 

where we can include the term R = in the sum because 
it vanishes. Using the orthogonality relation, we have 

^^.(k-k').R^^^^^,^^^^ (49) 

R K 

where the k are restricted to the first Brillouin zone and 
K are the reciprocal vectors of R satisfying 

K = 2fcA,m, (50) 

with m e Z3. Substituting Eq. in we obtain 
finally the expression [28| : 

P^^(k) = -nokfj^k^wCk) (51) 
-no ^ [{k^ + K^){k, + if,.)?I)(k + K) - K^K,w{K) 

where no is the number density of particles. In the grav- 
itational case, the integral (|45|l cannot be evaluated an- 
alytically. However, neglecting finite size effects, this in- 
tegral can be computed over the whole space and the 
periodic potential w(r) is approximated by the interac- 
tion pair potential w(r) = —Gm/r, so that 

w{\<i) ~ w(k) = / d?rv{v)e-'^ '' = (52) 

Using this it is straightforward to show (see App. the 
following simple result: 

3 

^w,^(k) = -nak^w{k) = inGpo. (53) 



In the context of the Coulomb lattice this is a well-known 
result, the so-called Kohn sum rule. In this case the 
quantity which appears on the r.h.s. of the sum, instead 
of AwGpQ, is —Wp — — 47re^no/TO where LUp is the plasma 
frequency. We will discuss further below the significance 
of this analogy. 

We can use these results and the above sum rule to 
compute — in a different way than in Eqs. (|2()|I - H21|I — 
the R = term of the dynamical matrix I?(R) (i.e. the 
term giving the force on a particle, at linear order in the 
relative displacements, when it alone is perturbed off the 
lattice). Using the Kohn sum rule (|53|l . the trace of the 
dynamical matrix is: 

tr[X>(R)] = AnGpo. (54) 

If the crystal has three equivalent orthogonal directions 
then the diagonal terms of the dynamical matrix will be 
equal. In the case of lattices with special symmetries 
(like the sc, bcc and fee) it is simple to show that when 
a single particle is displaced along the direction of an 
axis, the force acting on it is parallel to the direction 
of displacement . This implies that the non-diagonal 
terms of the dynamical matrix are zero. We can therefore 
conclude that 

25/..(0) = ^vrGpo^M- (55) 

2. The branches of the dispersion relation and the fluid 
limit 

We have noted that the spectrum of eigenvalues has a 
clear branch structure. To identify the different branches 
it is useful to consider the k ^ limit keeping the in- 
terparticle distance £ constant. We expect this to corre- 
spond to the fiuid limit: a plane wave fluctuation e*'' '' 
with k <C 1/^ has a variation scale much larger than the 
interparticle distance, and therefore does not "see" the 
particles. 

From Eq. 1)51(1 the limit for k — > is straightforward 
as the contribution of the sum on the r.h.s. vanishes in 
this limit 

lim Vf,iy(k) = -nofcpfc,yw(k). (56) 

k— *0 

Using the eigenvalue equation with Eqs. (|5T|l and 
()52|l . it follows that the solutions in the fluid limit are 

1. one longitudinal eigenvector polarized parallel to k 
with normalized eigenvalue £i(k ^ 0) = 1 and 



This can be explicitly shown e.g. using Eq. IA2I (taking the limit 
a — > and assuming that the sum over the replicas converges). 
We have assumed that the sum in Eq. 1511 is well defined — 
which is the case for the gravitational interaction — so that it is 
possible to take the limit before performing the sum. 



2. two transverse eigenvectors polarized in the plane 
transverse to k with normalized eigenvalues 
e2,3(k^0) = 0. 

As the spectrum of eigenvalues £„ (k) is exactly the same, 
up to an overall negative multiplicative constant, to that 
of the Coulomb lattice, we adapt the same terminology 
as in this context. The branch of eigenvalues whose asso- 
ciated eigenvectors converges to the longitudinal eigen- 
vector as k ^ is called the optical or longitudinal 
branch. The two other branches whose eigenvectors con- 
verge to the transverse eigenvectors are called the acous- 
tic branches. For finite k, the eigenvectors are not ex- 
actly parallel or perpendicular to k for all k but belong 
nevertheless to one of the three branches, which define 
three-dimensional hyper-surfaces in the four-dimensional 
space (w,k) space. 

The appearance of an optical branch in a monoatomic 
crystal is a characteristic feature of the 1/r interaction 
potential (at large r). In the case of a more rapidly de- 
caying potential at large scales, i.e., l/r^+" with a > 0, 
it becomes a third acoustic branch. In the case of a po- 
tential that decays slower at large r, i.e., a < 0, the 
optical branch diverges as k ^ 0. The physical inter- 
pretation of the optical branch is that it represents the 
coherent excitation of the whole lattice with respect to 
the background In a Coulomb crystal, the optical 

mode is produced by the lattice moving against this back- 
ground producing a "plasma oscillation" , at the plasma 
frequency ujp defined above. This mode is, as we have 
just seen, purely longitudinal, i.e., the perturbations are 
parallel to k, while the tranverse modes, i.e., the pertur- 
bations orthogonal to k have zero frequency. The reason 
for this behaviour of long wavelength density fluctuations 
can be easily understood. The density fluctuations are 
related, in this fluid limit, to the displacements through 
the continuity equation: 

(5/9 ~ V • u, (57) 

which implies in k space that 

- k • u. (58) 

Thus tranverse modes do not source density fluctuations, 
and therefore (by the Poisson equation) they do not pro- 
duce a force. In the case of gravity, instead of oscillating 
as in a plasma, the longitudinal mode may be amplified 
or attenuated (depending on the initial perturbation), in 
a way which is independent of fc. As we will discuss in de- 
tail below, this is just the well known linear amplification 
of density fluctuations in a self-gravitating fluid. 

3. Corrections to the fluid limit 

We have just seen that the fluid limit is obtained by 
taking the dynamical matrix as 

P(k) = ^k,k.. (59) 



— [1,0,01 




FIG. 3: Optical branch for different directions of k. The thick 
line is proportional to . 

We can estimate analytically the corrections to this limit 
for small k (i.e. for large wavelengths) by expanding the 
eigenvalues and eigenvectors of the full dynamical matrix 
about k = 0. We note that this corresponds to calcu- 
lating the difference, at large wavelengths, between the 
evolution of the perturbed lattice with a finite number of 
particles and that of the fluid limit. These are thus what 
are, in the context of cosmological simulations, "discrete- 
ness effects" introduced by the modelling of the fluid by 
such a system. We will discuss at length this application 
of this formalism in |Q . 

When expanding the dynamical matrix in Taylor series 
about the fluid limit k — > 0, it is simple to show that for 
1/r interactions this series is in even powers of k, because 
I?(R) is real and I?(k) analytic for k ^ (see jHls^). 
It is therefore possible to write the corrections to the 
eigenvalues of the optical mode as: 

u;f{k)~A7TGpo{l~h{k)k^), (60) 

where the expression for 5i(k) can be computed by di- 
agonalizing I?(k) expanded up to 0{k^). The leading 
correction to the two acoustic modes may be written 

w|(k) ~ 27rGpo62(k)fc^ (61a) 
w|(k) ~ 27TGpob3{k)k^. (61b) 

The Kohn sum rule implies that 6i(k) — (t'2(k) -I- 
63(k)) /2. In Fig. Owe show the optical branch, in various 
different chosen directions. The approximation with the 
leading term in the Taylor expansion is very good up to 
the Nyquist frequency. 

In Fig.2|we show how the anisotropy of the eigenvalues 
increases as the modulus of the wave vector increases (i.e. 
when we look at smaller spatial scales). We plot, for three 
ranges of values of the modulus of k, the value of the nor- 
malized eigenvalues as a function of the angle 9 between 
k and the axis that forms a minimal angle with it. As 9 
increases (i.e. as cos6' decreases with < 9 < 7r/2) there 
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(ii) 



FIG. 5: Schematic representation of (i) a mode collapsing 
faster than fluid limit and (ii) an oscillating mode. 



FIG. 4: Variation of the value of the eigenvalues for various 
ranges as a function of the cosine of the angle between k and 
the axes of the lattice which forms a minimal angle with it. 
We see that the effects of anisotropy are strongest for the 
short- wavelength modes, and decrease as we go towards the 
fluid limit. 

is a clear trend of decrease in the eigenvalue, in each of 
the three cases. The difference as a function of orienta- 
tion of the vector k is, however, much more marked for 
larger fc, i.e., at scales closer to the Nyquist frequency. 
This is not unexpected: the effects of anistropy (which is 
completely absent in the fluid hmit, in which the eigen- 
values are independent of the orientation k) are naturally 
strongest for the short wavelength modes. 



behaviour qualitatively different to that generically ex- 
pected based on the analysis of the fluid limit. Translated 
to the analagous Coulomb system, the result means that 
a sc lattice is, in this case, unstable (as there are grow- 
ing modes). While this result has not apparently been 
shown in the literature, it is not an unexpected result 
in this context. It has been established ,3^i42(i that for 
the (classical) Coulomb lattice that the ground state is 
the bcc lattice. It has a lower binding energy than the 
fee lattice, which in turn is a lower energy configuration 
than the sc lattice. Our result implies that the latter is 
not only a higher energy state, but that it is strictly un- 
stable. Indeed we note that the specific modes we have 
considered above describe a "sliding" of adjacent places 
in an sc lattice which deform it towards the lower energy 
configuration represented by the fee lattice. 



4- Oscillatory modes 



The spectrum of the sc lattice Fig. El includes some 
modes [e.g. for k = {kx,0,0)] with eigenvalues on the 
optical branch larger than the fluid limit. For exam- 
ple, this is the case for modes with initial displacement 
u(r,0) oc xexp{ikxx), shown in Fig. [Sj-f^ij. Adjacent 
planes collapse towards one another, faster than in the 
fluid limit. The Kohn sum rule Eq. H53|) states that the 
sum of the three eigenvalues a;^i(k) is equal to 47rG/3o. 
Therefore, the existence of modes collapsing faster than 
the fluid limit implies that there are other modes with 
negative eigenvalues 'j-'^(k), i.e., which oscillate. This 
is the case, e.g., of the mode with initial displacement 
u(r,0) ~ y exp{ikj;x), shown in the Fig.\^(ii). In this 
case, contiguous planes oscillate as indicated in the flg- 
ure. We will study these modes in greater detail using 
numerical simulation in Sect. IVT^ 

The existence of oscillating modes in a perturbed and 
cold purely self-gravitating system (i.e. without any 
additional interaction or velocity dispersion giving rise 
to a restoring pressure ^^) is an unexpected curiosity, a 



IV. GENERALIZATION TO AN EXPANDING 
UNIVERSE 

In the previous section, we have described the gravi- 
tational evolution of a perturbed lattice in a static Eu- 
clidean universe. In the cosmological context, density 
fluctuations are a perturbation around an homogeneous 
and isotropic Friedmann-Robertson- Walker (FRW) solu- 
tion of Einstein's field equations of general relativity. In 
cosmological N-body simulations, since the regions stud- 
ied are smaller than the Hubble radius and the velocities 
are non-relativistic, one considers the limit in which the 
equations of motion of the particles are strictly Newto- 
nian in physical coordinates r Q . These coordinates are 
related to the comoving coordinates x of the FRW solu- 
tion by 

r(t) - a{t)x{t), (62) 



If there is a non negligible velocity dispersion, it known that 



fluctuations at scales smaller than the Jeans length oscillate y|. 
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where a{t) is the scale factor describing the expansion of 
the universe. It satisfies the Fricdmann equation 



SttG 



-P- 



K 
9 ■ 



(63) 



where p is the mass density of the universe and k the cur- 
vature. In the unperturbed FRW model the particles are 
fixed in comoving coordinates, all deviation from these 
positions arising from perturbations to this model. For 
this reason it is very natural, and convenient, to work in 
comoving coordinates. We therefore start by transform- 
ing our previous Newtonian equations to these coordi- 
nates, the only further difference being that we perturb 
about a time-dependent solution describing an expanding 
FRW universe. 

Using Eq. H62|l the acceleration can be written 



r = ax 



2ax 



(64) 



The term ax can be expressed as the background con- 
tribution of the gravitational acceleration. For the spe- 
cific case of an Einstein de Sitter (EdS) Universe, i.e., 
a universe containing only matter without curvature 
[p{t) = p(j{a{t) / a{to)Y and n = 0], it is given by 

47r 

go = ax = — Gpox, (65) 
oa 

which has exactly the same form (for a = 1) as the con- 
tribution of the negative background of Eq. I|20() . We now 
write the position of a particle in comoving coordinates 
in terms of the displacement u about the lattice position 
as 



x(i) = R + u(R,i). 



(66) 



The vector R is now the position of the lattice sites in 
comoving coordinates (i.e. R does not depend on time) 
and u(R, t) is the displacement of the particle that was 
originally at R (in fluid theory, this is a Lagrangian co- 
ordinate^ see e.g. |2^). By using Eq. (|64|l . we can write 
Eq. H22|) in an expanding universe as 

ii(R,t) = -2-u(R,t) + ^ Vl?(R-R')u(R',t), (67) 
o?^ 

where we have implicitly included the background term 
(|65|l in the dynamical matrix. We emphasize that the 
dynamical matrix is identical to that in the static case: it 
depends only on the kind of lattice and on the interaction, 
but not on the dynamics of the background. Therefore all 
the analysis of this matrix performed in the preceeding 
section is valid also in this case. From Eq. H67|l . the mode 
equation H30(l generalizes simply to 



/„(k,i) + 2-/„(k,i) 
a 



(68) 



This is very similar to the equation of the evolution of 
a fluid in Lagrangian coordinates The difference is 
only in the factor u>n(k.) on the r.h.s., which in the fluid 
limit is replaced by AnGpo. 



A. Solution in an Einstein— De Sitter universe 

We derive now the solution of the mode equation H68() 
in the case of an EdS universe. The evolution of the scale 
factor is, from Eq. H63|l : 



ait) ^ ( |- 



2/3 



GnGpotQ — 1, 



(69) 



assuming that a(0) = 0. Then the mode coefficient equa- 
tion (IHHl) is 



/„(k,t) + -/„(k,t) 



3^2 



e„(k)/„(k,t). 



(70) 



where we have used again the adimensional quantity 
£n (k) defined in Eq. (|43|l . A set of independent solutions 
of ^7U^ which satisfies the IC are: 

-Q + (k) 



Un(k,t)^a{k) a+(k) 

F„(k,i) ^a(k)ta 
where 



to) 



a„ (k) 



'<(k) 



1 



an (k) -f a+(k) 



and 



a-(k) = - [v/l + 24e„(k)-l 



a+(k) = - + 24£„(k) + 1 
6 L 



(71a) 
(71b) 

(72) 

(73a) 
(73b) 



If £,i(k) > the solution presents a power-law am- 
plification mode and a power-law decaying mode. If 
— 1/24 < e„(k) < 0, there are two decaying modes. Fi- 
nally, if £,i(k) < —1/24, the solution is oscillatory and 
can be written as 



C/„(k, t) 





6 


a) 


cos 







7n(k) In 



1 



K(k, t) 
where 



67„(k) V^o 

to 



— sm 



7«(k)ln( - 



(74a) 



7n(k) 



— sm 



7n(k)ln - 
to 



7„(k) = -v/|24e„(k) + l|, 
6 



(74b) 



(75) 



i.e., the static oscillatory behavior of the static universe 
survives, but now the oscillation is periodic in the loga- 
rithm of time with decreasing amplitude. The evolution 
of the displacements is computed with Eq. (|34|l . The ef- 
fect of the expansion [through the "viscous" first term of 
the r.h.s. of Eq. (|67|l ] is to slow down the growing and de- 
caying mode of the non-expanding exponential solution 
into a power-law solution. 
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B. Fluid limit and Zeldovich approximation 

Let us calculate the fluid limit of the solution given by 
Eqs. 1(331), (123 and (ETJ- As explained in Sect. ITTll this 
corresponds to taking the limit k ^ at fixed £ of the 
dynamical matrix I'(k). In this case, as we have seen 
in Sect. IIIII one of the eigenvectors is parallel to k, with 
eigenvalue AirGpo, and the other two are normal to k 
with eigenvalue equal to zero. We have then: 

ei(k) = k, £i(k) = 1 — >ai ^ 2/3, ^ 1, (76a) 
e2(k) = k2±, e2(k) = ^ a+ 0, = 1/3, (76b) 
e3(k) = k3±, e3(k) = a+ = 0, = 1/3, (76c) 

where k2± and ksj. are orthogonal to k chosen so that 
k2_L • k3_L = 0. Using (|7()|l in (|7T|l . we get for the mode 
parallel to k: 



[/i(k,i)^C/||(i) = - 



3 / t 



2 V^o 
t 



2/3 



2/3 



, (77a) 
(77b) 



and for the modes perpendicular to k: 

C/2(k,i) = f/3(k,i) = C/i(i) = l, 



(78a) 



1/2 (k, t) = V3(k,i) = 1/^(0 = 3to 



1 



-1/3- 



(78b) 



The evolution operators (|35|l are then: 



P^^(k,t) = U\\{t)k^,K + (k2±);.(k2±)i. + (k3±)^(k3_L).., 

(79a) 

Q^,u{Kt) = V\\{t)kX+ (79b) 
+ VAt) [(k2±)M(k2±). + (k3±)M(k3±).] , 

[where we have used explicitly that U±{t) = 1]. Using 
Eq. (|34|l we write the evolution of the displacements in 
the fluid limit as: 

u(R,i) = Ui(R,io) + U||(R,to)C/||(i) (80) 
+v|| (R, to)Vi\ (t) + vx(R, to)Vx(t), 

where 

U||(R,io) - ^5]({i(k,io) •k)ke*-^, (81a) 

k 

uj.(R, to) = ^ 51 ("(1^' ^o) - ("(1^' ^o) ■ ^)k) e^*'-^, 

(81b) 



and analogously for the velocities v. Using the definition 
of peculiar gravitational acceleration g 



a 

—r = a 

a 



u 



2-u 

a 



(82) 



we can rewrite Eq. [with Eqs. (|73 and (O] as 
u(R,t) = Ui(R,to) 

+ g(R,to)io 



10 V*^ 



2/3 



V||(R,io)g^o 



1 



2/3 



3 /i_ 
5 (yi^^ 



-1/3 



+ Vi(R,io)3to 

where v is the peculiar velocity defined as 
v(x,t) 



a 

r r 

a 



(83) 



(84) 



This formula is precisely the one found in [29|. Note 
that in this reference, the Lagrangian coordinate X cor- 
responds to the position of the particle at t — Iq, i.e., 
X = R + u(R, to)- For asymptotically large times the 
solution is 



u(R,t).-to (- 



2/3 



-g(R,to)io +vj|(R,to) 



(85) 



This solution, using Eqs. (|82|) and H84|) . gives the follow- 
ing simple relation between the displacements and the 
peculiar velocity with the peculiar acceleration at any 
time: 



u(R,t) 



2 [to 
v(R,t) = g(R,io)i. 



4/3 



g(R,to)io: 



By imposing the IC 

u^(R,fo) = = v^(R,to), 

V||(R,io) = g(R,io)to = :5^U||(R,io), 

oto 



(86a) 
(86b) 

(87a) 
(87b) 



the evolution is given exactly at any time by Eqs. JHEJ, 
which is the well known Zeldovich approximation, in 
which the decaying mode is zero from the initial time. In 
cosmological N-body simulations of structure formation 
IC are canonically imposed 24., Al] using (|87|1 : given an 
initial power spectrum of density fluctuations a Gaussian 
realization of the gravitational potential is generated, and 
used to derive the initial gravitational field g(R, to) at the 
unperturbed particle positions. The particles are then 
displaced and given initial velocities as specified by H87|l . 
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COMPARISON WITH N-BODY 
SIMULATIONS 



In this section we compare predictions of the pertur- 
bative treatment we have presented in the previous sec- 
tions, which we will now refer to as particle linear theory 
(PLT), to what one obtains with full gravity (FG) cal- 
culated numerically in N-body simulations. The aim is 
to study the limits of validity of PLT as the relative dis- 
placements of the particles increase. We also compare 
with the evolution obtained using standard Lagrangian 
fluid linear theory (FLT). At the end of the section we 
also study further the oscillating modes which were iden- 
tified in the spectrum of eigenvalues discussed in Sect. lIIII 
Note that time will be expressed everywhere in this sec- 
tion in units of 1 / ^/AttGpq (i.e. the characteristic time 
scale in the fluid limit), unless otherwise stated. We also 
draw the reader's attention to the fact that the simula- 
tions we consider in this section are relatively small with 
respect to typical current standard cosmological simula- 
tions, but that this is not a relevant consideration here 
as PLT works for any finite number of particles. For our 
purposes here it is sufficient to have a large enough num- 
ber of particles to separate large scales (i.e. the box size) 
from small scales (i.e. the interparticle distance). 

The IC for our study are generated by perturbing 
in two different ways a sc lattice with N particles. We 
choose to restrict ourselves to simulations without space 
expansion and with zero initial velocities In this case 
it is straightforward to show, using the results of Sects. UTI 
and lllll that the displacements of the particles according 
to PLT are 



u(R,i) = ^^|exp(ik-R)x 

3 

A„(k) cosh ('V47rGpoe„(k) A e„(k) 



this becomes 



where en(k) and e„(k) are defined in Eqs. H28(l and (|43|l . 
while the A„(k) are determined by the IC. In the fiuid 
limit, following the discussion in Sect. IIIll after Eq. (|56|l . 



Except in the the study of oscillating modes in Sect. IV CI 
^'^ This choice is made for simplicity. Another equally simple case is 
that of the EdS universe, with velocities given by the Zeldovich 
approximation. The difference between the two cases is simply, 
as we have seen, in the time dependence of the modes. Indeed in 
the asymptotic fluid limit given by the Zeldovich approximation 
the two cases can be mapped onto one another by a simple trans- 
formation between time and scale factor: exp{v'47rGpot) i- In 
any case we will see that the criteria we deduce for validity of the 
PLT are expressed in a very simple way which one would expect 
to be essentially the same irrespective of the model. 



u(R,t) = l^|exp(zk-R)) 



^11 (k) cosh ( v47rGpo t k 



yl2i(k)k2_L + A3_L(k)k 



■3_L 



u(R,0) 



[cosh (V47rG>o t) - 1] 



u(R,0) . 

(89) 



This therefore corresponds to what we denote by FLT. 
Note that in the following, we consider this last equa- 
tion with the full initial acceleration u(R, 0), and not 
the linearized one [see Eq. (|13|l or (|16|l ]. since in the La- 
grangian fiuid approach, the force on a fiuid element is 
the full gravitational force [2^ . 

We consider two different kinds of initial displace- 
ments: spatially uncorrelated, and correlated. In the first 
case, each particle of the lattice is randomly displaced, 
with uniform probability, in a small cubic box centered 
on its lattice site. The power spectrum of density fiuctu- 
ation (PS) P(k, t) cx |(5p(k, t)p of the resulting particle 
distribution is proportional to fc^ at small k (421 143| . For 
the correlated case, the displacements are obtained from 
a set of Gaussian variables {a^(k), 6^(k)}j^ ^ : the /xth 
component of the displacement of a particle at the lat- 
tice point R is then 



(R) = 1 ^exp(zk • R) K(k) + tb^{k)] 



(90) 



The random variables a^(k) and bf_i(k) have average 
and variance cr^(k) cx k~'^ and are statistically indepen- 
dent of one another This gives rise, to a good approx- 
imation, to a distribution with PS of the density field 
proportional to fc~^ at small k 

We have performed our FG N-body simulations using 
the publicly available treecode Gadget j4^] In ta- 
ble ^ are summarized the parameters characterizing the 
two different IC and simulations. The number of parti- 
cles in the simulations is 16^ and 32^, respectively. Note 
that the softening lengths e chosen are in both cases 
much smaller than the initial inter-particle distance. This 
means that this modification of the gravitational force 
can be neglected in our calculation of PLT since it has 



Note that for it^{R) to be real, the following conditions are re- 
quired a^(-k) = a^(k) and 6^(-k) = -f)^(k). 
This can be seen directly from Eq. 1581 . which is valid in the 
continuous limit, i.e., in the limit of small perturbations and 
when the effects of the discretisation introduced by the lattice 
are neglec ted . For a detailed analysis of the corrections to this 
result see EE Es, 44]. 
^"^ Version 1.2. See http:/ /www.mpa-garching.mpg.de/gadget/] 
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Name 




P{k) 




e/ dNN 


Cd(^,0) 


k 


UC16 


16 




0.97 ±0.01 


0.029 


0.002 


5.0 


C032 


32 




0.983 ± 0.009 


0.057 


0.0006 


5.0 



TABLE 1: Summary of the two different distributions used as 
IC for the comparison between hnear theory and true gravity 
(UC stands for uncorrelated while CO is for correlated). The 
column " dNN /£ " gives the average distance between nearest 
neighbors and its standard deviation, both normalized by £. 
The next column is the softening length e, normalized by dNN, 
used in the FG simulation. The fifth column gives the initial 
value of the function ("_d at the interparticle distance, defined 
in Eq. 19111 . The time in the last column is defined later in 
the same equation: it gives the time scale at which PLT is 
expected to break down. 



already broken down when nearest neighbor particles are 
so close For both sets of IC, we have considered three 
different evolutions: (i) according to FG, (ii) according 
to PLT and (iii) according to FLT. Snapshots are shown 
in Figs.Eland[71 In these projected representations of the 
configurations at different times, we have considered in 
the latter case (C032) a subvolume of the simulation box 
containing initially 16^ particles in order to facilitate the 
visual comparison of the two cases. From these figures 
one can already see how well, to a first approximation, 
PLT works in predicting the evolution of the two sys- 
tems considered. The improvement given by PLT over 
FLT as an approximation to the FG evolution is visually 
very clearly manifest in the case of UC16. In the case 
of C032 the difference between the PLT and FLT is less 
visually evident, but one can still discern that the former 
does distinctly better in tracing the non-linear structures 
forming in the FG simulation. The reason why the differ- 
ences are more pronounced in the case of UC16 is simply 
that in the correponding IC there is, compared to C032, 
relatively much more power at small scales than at large 
scales. As we have seen in Sect. lIIll it is at these smaller 
scales that the two approximations differ. 

The amplitudes of the initial displacements in both 
simulations have been chosen sufficiently small so that 
linear theory remains valid during a non-negligible time 
(in units of l/V47rGpo). Actually both Eqs. (jHEl and lIH?^ 
predict the same evolution when t is small: u(R, t) — 
(t^/2)u(R, 0). In order to distinguish predictions of PLT 
from those of FLT, the linear regime must last sufficiently 
long, i.e., the initial relative displacements of particles 
must be sufficiently small. We will see below that, for 
our chosen initial amplitudes, PLT breaks down approx- 
imately at the same time in the two systems. This is 
purely coincidental. 



The effect of smoothing, which can be easily implemented in 
the PLT, will be discussed in a forthcoming paper 23j on the 
quantification of discreteness effects. 



Let us first consider how to characterize quantitatively 
the regime of validity of PLT. As discussed in Sect. Ill Bl 
we anticipate that the relevant quantity will be the rel- 
ative displacement of particles compared to the initial 
distance separating them. This motivates the study of 
the following statistical quantity: 

(|u(R',t)-u(R- + R,t)p)^, 

C,D[ti.,t) = p— ^ 

l^^l (91) 

- ^(eD(o,i)-eD(R,t)) 

where ^d(R, t) = (u(R', t) ■ u(R' + R, t))^, is the trace 
of the correlation tensor of the displacements which is re- 
lated to the PS of displacements PD(k,t) = |u(k,t)p/iV 

by 19 

(R, = ^ E ^) • ^) • (92) 

k 

Given the convergence criterion Eq. (|15|l . we expect PLT 
to break down when Cd (R, t) 1 for some R. One would 
expect that this condition will first be attained for |R| = 
£, i.e., when a significant number of nearest neighbors 
have come close to one another. 

Figures |H1 and El show the evolution of the function 
Cd(R, t) averaged over the first six nearest neighbors, i.e., 
for the six vectors R of norm £ 2°, for IC UC16 and C032 
respectively. In both cases, the FG evolution is compared 
with the predictions of PLT and FLT. The horizontal line 
is at 1 and allows one to determine the characteristic time 
t,^ such that 

CD(Atc) = ^ E CD(R,ic) = l> (93) 
\R\=e 

around which we expect PLT to break down, as discussed 
above. This time — reported in tab. HI — is of the order 



We note that for the case of the spectrum of displacements 
in C032, P£)(k, 0) oc |k|~*, this quantity is actually not well 
defined, in the sense that it diverges in the thermodynamic 
limit (i.e. L ^ oo at fixed number density no = N/L^). 
Indeed for a generic PS of displacements PoO^) <x |k|", we 
have, taking the continuous limit of Eq. I92i . that ^d(R-) = 
2Ty'i\i R Jo" ^"^^ sin(A;i?) dk which is infrared divergent for n < 
— 3 (kc being an ultraviolet cut-off given by the Ny quist frequency 
in this context). However, it can be shown Il5l that the func- 
tion fj5(R, t) is well defined in the same limit for n > —5. This 
limit in fact just coincides with the condition that the associ- 
ated density fiuctuations have finite variance, since this requires 
that limk^Q k^ {k^ PDCk.)) = 0. These divergences in ^d{R-) are 
not a problem provided (R, t) is well defined: the particles 
can move an infinite distance from their lattice positions, but 
what is important for the validity of the approximation used is 
how much their relative displacements changes compared to their 
separation. 

Three of them are actually enough since (R, t) is symmetric 
in its first argument. 
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FIG. 6: Evolution from IC UC16 (projection on the plane z = 0). From left to right, times 4.5, 5 and 5.5 and from top to 
bottom, FG, PLT and FLT. At f = 5, PLT breaks down according to the criterion based on the function (u{i,t) defined in 
Eq. H93|l (see the time in tab.|l|. 



of 5.0 for both IC We see from these figures that this 
expectation turns out to be correct: the PLT evolution 
from both sets of IC traces very accurately the FG evo- 
lution of the N-body simulation until a time very close to 
tj. For UC16 (Fig.lSJl the deviation between the PLT and 
FG evolutions becomes significant for a time just slightly 
earlier, while for C032 (Fig.O this time is slightly later. 
Considering the curves for FLT, we see that it does less 



As mentioned above, the fact that this time is very close in both 
simulations is purely coincidental. 



well than PLT in both cases, the improvement given by 
PLT over FLT being much more marked in the case of 
UC16. This confirms quantitatively the visual impres- 
sion of Figs. El and 13 for the same reason we gave above. 
Even though (Y){£,t) is a real space quantity character- 
ising correlation around the inter-particle distance i, it 
is an integral in k space which picks up a contribution 
from longer wave-length modes which are well described 
(as discussed in Sects. ^ and IIIII and further below) by 
FLT. For the much more infra-red dominated spectrum 
of displacements of C032 this contribution is much more 
significant and so FLT naturally is a better approxima- 
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FIG. 7: Evolution from IC C032 (projection on the plane z — 0). From left to right, times 4.5, 5 and 5.5 and from top to 
bottom, FG, PLT and FLT. At f = 5, PLT breaks down according to the criterion based on the function (-D{£,t) defined in 
Eq. H9.'i^ (see the time in tab.QJ. Note that this is a cubic subset of 1/8 of the system to have approximatively 16'^ particles, 
as for UC16 in Fig. H 



tion for this quantity than in the case of UC16. 

We next give an alternative, perhaps intuititively more 
direct, way of quantifying the regime of validity of 
PLT. We consider the nearest neighbor (NN) distribution 
u}{r,t): at a given time t, this function gives the proba- 
bility density for a particle to have its nearest neighbor 
at the distance r (see e.g. ^3)- Figs. ITUl and ITTI are 
shown, for the FG evolution of UC16 and C032, the cu- 
mulative distributions derived from ti;(r, t), 



(94) 



T(r,i) = / uj{s,t)ds 



These represent the probability that a given particle has 
its NN within a distance r. It allows one to determine 
quite accurately the time at which "shell crossing" oc- 



curs . We see that between i « 4 and t 



t(^, the 



^■^ We adopt here loosely the terminology used in fluid theory to 
refer to the time when particles fall on top of one another: if one 
considers the particles in the simulation as the centres of fluid 
elements, this corresponds to what is called "shell crossing", at 
which point the linearised Lagrangian fluid theory (i.e. FLT) 
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FIG. 8: Evolution of the function (,o{£, t), defined in Eq. (CT . 
according to EG, PLT and FLT in the UC16 simulation. The 
horizontal line is 1 and gives the characteristic time ti^ around 
which we expect PLT to break down (as an approximation to 
EG). 
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t 



EIG. 9: Evolution of the function (-D{£,t) according to EG, 
LPT and ELT in the C032 simulation. The horizontal line is 
1 and gives the characteristic time around which we expect 
PLT to break down (as an approximation to EG). 



behavior of the cumulative distributions changes in an 
important way: for t < A, almost no particle has its NN 
closer than a distance r « 0.3, while at time i w 5 = i^, 
50% of the particles have their NN at a distance smaller 
than this distance. Thus at i w 4, the first shell crossings 
occur and at i « 5, approximately half of the particles 
have already had their own shell crossing or are very close 
to it. 



breaks down as the density diverges. 
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EIG. 10: Evolution of the cumulative NN distribution T(r, t) 
in the UC16 simulation. The times are indicated in the legend. 
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EIG. 11: Evolution of the NN distribution T{r,t) in the C032 
simulation. The times are indicated in the legend. 



A. Comparison of relative displacements 

It is interesting to study also the accuracy of PLT in 
tracing the FG evolution of the trajectories of individ- 
ual particles (rather than only averaged quantities). We 
consider now the relative displacement of a particle with 
respect to its NN, that is u(R, t) — u(R', t) where R' and 
R are separated by a vector of elementary size i. To do 
so we have chosen randomly a particle in each simula- 
tion and selected among its six initial NN the one which 
ends up closest to it at the time at which PLT breaks 
down in the corresponding FG simulation (see Figs. El 
and [T^ . In Figs. and El are shown the result for 
UC16 and C032 respectively. These two figures allow 
one to see that for both IC, PLT describes very well the 
evolution of the relative displacement of the two particles 
considered up to a time of the order of 4.5 = 0.9t^. In 
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FIG. 12: Evolution of the distances d{t) (normalized to £) 
between a randomly chosen particle in UC16 and its six NN 
in the FG simulation. The arrow shows the NN which is used 
for the study of the relative displacement (Fig. 114^ . 



2 



1.5 - 




I ' ' ' ' ' ' ' ' ' ' ' ' 1 

1 2 3 4 5 6 

t 

FIG. 13: Evolution of the distances d{t) (normalized to £) 
between a randomly chosen particle in C032 and its six NN 
in the FG simulation. The arrow shows the NN which is used 
for the study of the relative displacement fFig. I15t . 



both systems, this time is slightly smaller than the time 
at which the particle selected is the closest to its NN: 
from Figs. 1121 and 1131 we may estimate that this time is 
approximately 5.0 in the system UC16 and 5.5 in C032. 

We also note that PLT does again also in this case 
better than FLT in describing the relative displacements. 
Indeed FLT already breaks down at t 3 in UC16 and 
t « 3.5 in C032. 



B. Evolution of modes 

In this subsection we consider the comparison of PLT 
with FG and FLT in reciprocal space. We first consider 



2.5 r- 

2 - 




-0.5 



1 2 3 4 5 



FIG. 14: Evolution according to FG (thick line), PLT 
(medium line) and FLT (thin line) of the absolute value of 
each coordinate (denoted by a different type line) of the rel- 
ative displacement u(R, t) — u(R',t) of a randomly chosen 
particle and its NN in UC16. The particle is the same as the 
one chosen for Fig. 1121 and the NN corresponds to the one in- 
dicated by an arrow in this last figure. Note that, for clarity, 
we have shifted two of the coordinates by ±0.4. 
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FIG. 15: Evolution according to FG (thick line), PLT 
(medium line) and FLT (thin line) of the absolute value of 
each coordinate (denoted by a different type line) of the rel- 
ative displacement u(R, t) — u(R',t) of a randomly chosen 
particle and its NN in C032. The particle is the same as the 
one chosen for Fig. 1131 and the NN corresponds to the one in- 
dicated by an arrow in this last figure. Note that, for clarity, 
we have shifted two of the coordinates by ±0.3. 



the evolution of the PS of displacements PuO^, t) for a few 
specific vectors k. Then we study the evolution of this 
quantity, but now averaged on vectors k of similar modu- 
lus. Besides verifying the conclusions drawn in Sect. ITm 
this numerical study allows us to assess the validity of 
PLT for different wavenumbers as a function of time. 
Following the results and discussion in Sect. IIIII we 
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expect firstly to see the PLT evolution to differ less and 
less from FLT as we go to longer wavelength modes 
(|k| <C 1/^), since the eigenvalues and eigenvectors of the 
modes approach those of the fluid limit in this case. For 
such modes we expect that PLT and FLT should both 
follow the FG evolution accurately up to at least t ^ t(^. 
For short wavelength modes (|k| ~ l/£)we expect PLT to 
be significantly more accurate than FLT. As for the time 
of breakdown of PLT for a given mode, we would expect 
that long wavelength mode evolution should be described 
accurately by PLT for a time longer than that of short 
wavelength modes: the breakdown of the approximation 
in the sense we have characterised it above, i.e., in terms 
of the approach of NN particles, would not be expected 
to affect significantly the evolution of longer wavelength 
modes. We will discuss this point further below. 

Figures ITHI and [T7I show the evolution of PD(k, <), nor- 
malized to its initial value (i.e. at t — 0), for two chosen 
vectors k with very different lengths and different ori- 
entations with respect to the lattice, for IC UC16 and 
C032. In both cases, PLT follows very accurately the 
FG evolution up to a certain time for both long and small 
wavelength modes. The time up to which the agree- 
ment is good depends, as has been anticipated, on the 
wavelength. It breaks down first for the large k (small 
wavelength) mode, and significantly later for the small k 
(large wavelength) mode. For both IC the effects of non- 
linearity become significant for the large k mode chosen 
(one of the modes of largest modulus in each case), at 
t « 3.5, i.e., slightly before the time of the first shell 
crossings as determined above from Figs. 11111 and ITTI In 
this case we see also clearly that FLT is a poor approx- 
imation to the evolution, corresponding to an evolution 
with an exponent which is significantly too large. For 
the long wavelength mode (one of the modes of smallest 
modulus in each case), on the other hand, we see that 
FLT and PLT, as expected predict the same evolution 
(at the level of precision allowed by the figures). And in 
this case we see that they both trace the FG evolution 
very accurately for a very significantly longer time, up to 
t « 6 in the case of C032. 

We now consider the PS of displacements PdO^, t), av- 
eraged in spherical shells, which we denote simply by 
Pu{k,t): 

Pj,(k,t) = — y PD(k,t) (95) 

fc-dfc<|k|<fc+dfe 

where is the number of vectors k in the spherical shell 
[A; — dfc, k -|- dfc[. This analysis of Poik, t) allows one to 
observe in further detail how PLT fails in describing FG 
evolution after some time and at different scales. More- 
over, it shows how it is more accurate than FLT in ap- 
proximating FG. 

Figure^]shows the evolution of the averaged PS of dis- 
placements according to FG (thick lines), PLT (medium 
lines) and FLT (thin fines) with IC UG16. The first differ- 
ence between PLT and FG appears at large k&ttK, 4.5, 
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FIG. 16: Evolution of PDu(k, t) normalized to its initial value 
for two chosen vectors k according to FG (thick lines), PLT 
(medium lines) and FLT (thin lines) with IC UC16. The 
vectors k are 27r(l, 0, 0) /L and 27r(7, 7, 7) /L. 
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FIG. 17: Plot similar to the one in Fig. [HI but for IC C032 
and vectors k = 27r(l,0,0)/L and k = 27r(15, 15, 15)/L. 



slightly before i<^. This difference propagates to smaller k 
at later times. FLT is already discernibly different from 
FG at i w 1.8 at large k and the difference propagates to 
smaller k at later times. 

Figure [T^is similar to Fig. [TBI but concerns IC G032. 
Note that, since PD(fc,0) cx; fc^"*, it is Pd(^,0^^ which 
is shown Similar conclusions can be drawn to those 
in the uncorrelated case. The differences between PLT 



A slight excess of power over the expected PS is clearly visible 
around the Nyquist frequency. This is a small aliasing effect due 
to the fact that we have included some k in our sum outside the 
first Brillouin zone: we have summed in Eq. 1901 over the modes 
(L/27r)k e [-30,30]^ rather than (L/27r)k e [-16, ISf. This 
has no bearing on the conclusions drawn here. 
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FIG. 18: Evolution of the averaged PS of displacements 
Pu{k,t), Eq. according to FG (thick lines), PLT 

(medium lines) and FLT (thin lines) from IC UC16. The 
times are indicated at the left of the curves. The values of k 
are given in unit of the Nyquist frequency /cn- 
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FIG. 19: Evolution of the averaged PS of displacements 
Po{k,t), Eq. according to FG (thick lines), PLT 

(medium lines) and FLT (thin lines) from IC C032. The 
times are indicated at the left of the curves. The values of 
k are given in unit of /cn- Note that it is actually Pi){k,t)k'^ 
which is represented. 



and FG become visible at i w 4.5 and propagate at later 
times. Actually at t « 6.5, there is no longer good agree- 
ment at any k. FLT starts to deviate discernibly from 
FG at i « 1.8. 

Figs. 1181 and 1191 thus confirm clearly what was already 
observed in Figs. ^1 and El the breakdown of PLT 
starts at the largest k and propagates progressively in 
time to smaller k. Further PLT (and FLT) remain a 
good approximation to the evolution at the smaller k 
in our simulations at times significantly longer than the 
time ti^ which we used to characterize the global valid- 



ity of PLT. The reason is simple^'*, and it is the same 
one which explains, e.g., why linear fluid theory success- 
fully describes the evolution of small perturbations to a 
self-gravitating system even when there are strong non- 
linearities at smaller scales: re-arranging matter in any 
way, subject only to the constraint that matter and mo- 
mentum are conserved, up to some real space finite scale, 
£c say, can produce, at most, fluctuations at small k (i.e. 
k l/£c) with a PS of density fluctuations cx /c"*. The 
perturbative approximation to full gravity represented 
by PLT breaks down globally, as we have seen, when 
NN start to approach one another. When this happens, 
however, the full gravitational force can still be approx- 
imated for a longer time by the full gravitational force 
on each particle due its NN particle, plus the force from 
all other particles still linearised, as in PLT, in the rel- 
ative displacements from their starting positions. This 
means that FG can continue to be approximated by PLT 
plus an effective interaction of finite range. This latter 
interaction can produce at most a term in the PS in fc*, 
which will always be dominated by the simple amplifica- 
tion given by PLT of the initial PS, which in the cases we 
have considered has a small k behaviour cx fc" with n < 4 
(n = 2 and n = — 2 for UC16 and C032 respectively). 
An interesting question is whether, at these longer times, 
PLT is a good approximation in so far as it agrees well 
with FLT, or whether it can actually continue to trace 
the FG evolution better than FLT. From the numerical 
results we have given here it is not possible to distinguish 
FLT from PLT in the corresponding regime. Larger sim- 
ulations would be required to answer this question. 



C. Oscillating modes 



We have noted in Sect. Illll that the spectrum of eigen- 
values of the sc lattice contains some negative eigenval- 
ues, which give oscillating modes. These modes are not 
of practical importance in typical cosmological simula- 
tions since they are relatively few, and in a short time of 
little importance compared to the unstable modes (which 
also have considerably larger exponents ^^). It is never- 
theless interesting to study them briefly since they are a 
peculiarity of the discrete system. Indeed, as discussed 
in Sect. IIIII such oscillating modes do not exist in fluid 
theory without initial velocities. 

To study these modes we consider a sc lattice with 



^'^ The argument given is, in the context of cosmology, attributed 
to Z eldovich. For an extensive discussion see 1| §27-28, and also 

m. 

Further, the IC of cosmological simulations, because of the Zel- 
dovich approximation Eqs. I85I - IS71 . are purely longitudinal, 
while the oscillating modes are on the accoustic branches which 
are close to purely transversal for most k. 
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iVl/3 



even and the following initial displacement: 



u(R, t) = 



+(5y if Rx/t is even, 
— (5y if Rx/i is odd, 



(96) 



where (5 is a constant. This corresponds to displacing 
each plane of particles with constant R^ in the direction 
+y and — y alternately (see Fig. |3Jl. The only excited 
mode is k = (0, — 1/£, 0). According to PLT, the eigen- 
value associated to this displacement is —0.156 • 47rG'/9o/3 
(cf. Sect. lIII|l and the particles therefore should oscillate. 
If no initial velocity is considered, then the motion of a 
particle should be, as predicted by PLT 



u(R, t) = ±(5 cos (cjQ t) y, 



(97) 



with a frequency loq — ^0.156 • 47rG/3o/3. 

To observe these oscillations numerically, we have writ- 
ten a special code to integrate the equations of motion 
with FG, rather than using the same code (Gadget) as 
in the preceeding subsections. The reason we do this is 
that it is very difficult numerically to observe them with 
such a code. Indeed the simulation of oscillating modes 
would be an interesting and challenging test for the pre- 
cision of gravitational N-body codes. The primary diffi- 
culty is that the negative eigenvalues are small compared 
to most of positive ones (cf. Fig.|21in Sect. lIIip . so that an 
even much smaller amplitude perturbation in any of these 
modes can grow as an instability on a time scale much 
smaller than the period of oscillation. Such small pertur- 
bations are clearly created by numerical imprecision, but 
also, by the intrinsic non-linearities in the FG, i.e., there 
is a coupling of modes which is neglected in PLT. Thus 
one must work not only with great numerical precision 
but also at extremely low amplitude to avoid "contam- 
ination" by other modes through the non-linearities on 
the relevant (long) time-scale. 

The code that we have used has been built specifically 
for the particular initial displacements H96|l taking into 
account the following considerations. From the symme- 
tries of the configuration, it follows that the full gravi- 
tational force on a particle is only along the y axis, and 
modulo a change of sign, all the particles move the same 
way. Moreover, since the distribution of particles can 
be seen as two perfect rectangular sub-lattices with lat- 
tice spacing (. in the y and z directions and 2€ in the 
X direction the force on a given particle is only due 



Some eigenvectors of the oscillating modes have in this case the 
particularly simple form we will treat. The fact that these modes 
differ so significantly for the case that is odd or even illus- 

trates again that their presence is entirely associated with the 
discreteness of the system. 
^"^ One lattice is constituted by the particles which are initially dis- 
placed by an amount +i5 and the other one by the particles dis- 
placed by —(5. 
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FIG. 20: Oscillations of a particle with initial displace- 
ment 1961 1 and 5 — 0.004£ according to FG (solid line) and 
PLT (dashed line). The curves are actually indistinguishable 
on the plot. The time is in units of l/y/^nGpo and the dis- 
placement is normalized to the value of S. Details of the FG 
simulation are given in the text. 



to the particles in the sub-lattice which does not con- 
tain this particle In order to evolve the system, it is 
therefore sufhcient to integrate the equation of motion 
of a single particle in one dimension with a force com- 
ing from half of the particles of the system (the force is 
calculated by using the Ewald summation formula) . The 
method used for the integration is the embedded Runge- 
Kutta-Fehlherg (4, 5) method, implemented in the GNU 
scientific library and more precise than the standard 
leap-frog method used in cosmological N-body simula- 
tions. One can also note that, due to the periodicity of 
the system'^", the number of particles is not important 
as long as N^/^ is even: N^/^ = 2 is actually enough. 

Figure 1201 shows the oscillations along the y axis of a 
particle according to FG (obtained by using the code we 
have just described) and the prediction of PLT. The value 
of 5 used is 0.004 ^, for which we find a perfect agreement 
between FG and PLT. It is interesting to note that for 
larger value of 5, the frequencies of the oscillations mea- 
sured in the FG simulations decrease and the functional 
behavior of the oscillations is less and less close to ex- 
actly sinusoidal as can be seen in Fig. [2] This trend 
towards a decreasing frequency [rather than the constant 
frequency ojq as in Eq. H97|l ] as the amplitude increases 
has a simple explanation: when (5 = ^/4 the full force 
exactly vanishes as the perturbed configuration is in this 
case again a perfect sc lattice. This is also true when 



^® The force from the particles in the same sub-lattice is zero by 
symmetry. 

http:/ /www. gnu.org/software/gsl/. 

The systems which we consider are always periodic at the level of 
the box but here the distribution inside the box is itself periodic. 
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FIG. 21: Oscillations of a particle with initial displace- 
ment II96II and & = 0.248^ according to FG (solid line) and 
PLT (dashed line) . The last curve ( "Best fit" , dashed-dotted 
line) represents a sine with frequency ^0.0116 • 47rG'po/3. It 
allows one to see that the frequency of the oscillations is 
smaller than the one predicted by PLT and that the func- 
tional behavior is slightly different from an exact sine. The 
time is in units of 1/ V47rGpo and the displacement is normal- 
ized to the value of &. Details of the FG simulation are given 
in the text. 



(5 — 1/2, but in that case the resulting distribution is just 
the initial sc lattice. It follows that if 5/€ e]l/4, l/2[, one 
observes the same types of oscillation as for 5/^ e]0, l/4[ 
due to the invariance of the system under transformation 
of the type 5 jl/2 ± 5 with j an integer. 



VI. CONCLUSIONS 

In this paper we have described in detail a new pertur- 
bative treatment which describes the evolution of N self- 
gravitating particles of equal mass initially perturbed off 
a perfect lattice, subject to periodic boundary conditions, 
both in a static spacetime and in a cosmological (expand- 
ing) background. We have reported specifically the spec- 
trum of eigenvectors and eigenvalues for the modes of 
the displacement field on a simple cubic lattice, which is 
the case of relevance in cosmology. While the fluid limit 
[N oo) is recovered for long wavelength modes, the full 
spectrum for the finite N system contains both modes 
with negative eigenvalues, corresponding to oscillations, 
and modes with exponents greater than in the fluid limit. 
Further the eigenvalues depend explicitly not only on the 
modulus of the wave- vector k, but also on its orientation 
with respect to the axes of the lattice. The breaking 
of rotational invariance in the lattice is thus imprinted 
in the evolution of the system. We have shown, by de- 
tailed comparison with numerical simulations, that the 
linear order of the scheme has approximately the same 
range of validity as the corresponding (linear) order of 



the fluid theory, up to when particles come very close 
to one another (i.e. up to "shell-crossing" in fluid lan- 
guage) . However it traces the real evolution with greater 
accuracy than its fluid counterpart. 

In the context of cosmological simulations this means 
that our method provides a precise tool for quantifying 
fully, up to shell crossing, the effects of discreteness in 
these simulations: these effects are nothing other than 
the difference between the finite N evolution and the fluid 
limit. In a forthcoming paper |23| we will explore more 
extensively this application, giving precise quantitative 
measures of these effects adapted for use in "correcting" 
such N-body simulations. We conclude this paper by 
commenting on a few possible developments of the per- 
turbative theory we have described here. 

Firstly the method can evidently also be extended to 
higher order, just as has been done in the analogous treat- 
ment of condensed matter system (see e.g. ^3;I13)- It 
is straighforward to generalize Eq. IjlGII to an expansion 
to any order. The /i component of the force reads: 



= E E (98) 

11=0 R^R. 

X (R) - u^, (R')] . . . [u^„ (R) - (R')] ■ 



The tensor C/^"2i...iy„ is a function only of the interaction 
potential v{r) and is given by: 



a("+iV>(r) 



(99) 



r=R 



The analysis can be followed through at any order in 
analogy to linear order. Transforming to reciprocal space 
the problem simplifies, but there is now the added com- 
plexity of a coupling of modes. A first calculation of 
interest would be to map this description, in the fiuid 
limit, onto the corresponding one at the same order in 
the fluid Lagrangian theory. This latter treatment has 
been explored extensively in the cosmological literature 
and compared in detail with numerical simulations (see 
e.g. and references therein). The study of the cor- 
rections to this limit should allow one to get some insight 
on the interplay of non-linearity (which is of course also 
a feature of the fluid model) and the effects of discrete- 
ness. Non-linearity in gravity involves the transfer of 
power from larger to smaller scales, an effect which is 
often qualitatively argued to make discreteness effects of 
less consequence. One might hope to see such a mecha- 
nism at play, if indeed it is there. 

The approach presented here may also prove useful in 
providing insight about the nature of existing approxi- 
mations to self-gravitating systems which go beyond the 
simple fluid limit, as it provides an "exact" evolution of a 
self-gravitating flnite N-body system in a certain range. 
For example approximations have been developed to self- 
gravitating systems involving pressure terms associated 
to velocity dispersion (see e.g. El El El HI IH HI ) . 
In principle these terms can be calculated exactly using 
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our perturbative scheme and the improvement (if this 
is the case) they allow to the approximation of the full 
evolution better understood. 



Another direction in which this treatment can be gen- 
eralised is to the consideration of other initial configura- 
tions. We have analysed here almost exclusively the case 
of perturbations from a sc lattice, as this is the kind of lat- 
tice used in cosmological simulations. Our treatment can 
easily be generalised to other lattices, and a comparative 
study of the discreteness effects should be straightfor- 
ward. Without doing any calculation however one simple 
and interesting result can be given. It is known that both 
the bcc and fee lattice are stable (or at least meta-stable) 
configurations for the Coulomb lattice. This means that 
there are no unstable modes. For gravity this implies 
that there are only unstable modes. There are therefore 
no oscillating modes, and thus, by the Kohn sum rule, 
no modes with exponent greater than in the fluid limit. 
Consequently either of these lattices would appear to be 
better lattices to use in N-body simulations, in which one 
wishes evidently to approach as closely as possible the 
fluid evolution. The bcc lattice would appear to be the 
more interesting of the two, as it is known to be the 
most densely packed lattice. It is also a more isotropic 
configuration than either the sc or fee lattice. 



Another case of interest to analyse would be that of 
"glassy" configurations, which are often used as an al- 
ternative to the sc lattice in cosmological simulations 
0, li^ . These configurations are generated by sim- 
ulating a set of point particles evolving under negative 
gravity (i.e. Coulomb forces), with an appropriate damp- 
ing term. By doing so one arrives at a configuration 
in which the forces are very small, but which is more 
isotropic than the lattice (in which the forces are exactly 
zero). In the approximation that the initial forces are 
negligible one could, in principle, carry out the same kind 
of analysis as given here. The only difference is that the 
37V eigenmodes of the displacement field will not be plane 
waves, which greatly complicates the analysis compared 
to the case of the lattice. Numerically however such a 
solution should be feasible (for any specified glassy con- 
figuration), and it would be necessary if this method is 
to be used to give a precise quantification of discrete- 
ness effects from these IC as we can now give, using the 
analysis presented here, for the case of a lattice. We do 
not expect the results, however, to be very different (ei- 
ther qualitatively or quantitively): the effects described 
here are essentially sampling effects which depend on the 
sampling scale (the lattice spacing £ above) and not on 
the precise nature of the sampling. The particular mani- 
festation of anisotropy which we have observed in the sc 
lattice will necessarily be quite different, and likely less 
pronounced, in the glassy case, but the average slowing 
down of growth at smaller scales would be expected to 
be very similar in magnitude. 
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APPENDIX A: EWALD SUM OF THE 
DYNAMICAL MATRIX V{R) 

The Ewald sum for the dynamical matrix is given from 
(|T7|l using the Ewald sum for the potential : 



(Al) 



with 



(R- nL)^(R- nL) 



|R-nL|' 



4Q! 



X ^exp(-Q;^|R- nLp) 



|R-nLP 



-3 



(R-n£)^(R-nL), 
|R- nL|5 



(A2) 



2a 



erfc(a|R - nLl) + —= exp(-a2|R - nLP)|R - nL\ 



and 



47rGm 



^ in; 



|k|2\ 

exp [ —-^ ) cos (k • R) fc^fc^. 



4a2 ; 



The R = term is 



I?(R = 0) = - ^ I?(R). 



(A3) 
(A4) 



Note that, by symmetry, only the first term of the r.h.s. 
of HA2p and Eq. (|A3|) contribute in the sum of Eq. IIA4p . 
In the case of pure gravity the result of the sum IIA4II is 
given by Eq. ((^OJ) . 



APPENDIX B: KOHN SUM RULE 

We derive here the Kohn sum rule (|53|l . Multiply- 
ing Eq. H51|) by (e„(k))p(e„(k))iy and summing over 
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n, fi and v we obtain, with Eq. H28|l : 



+ ^?I;(k + K) [(k + K) •e„(k)]' 
-^^(K)[K.e„(k)]4. (Bl) 



Using the orthogonahty relation 

3 



^(e„(k))^(e„(k)).=<5, 



(B2) 



n=l 

we get finaUy ^ 

3 



(B3) 



-no ^ (|k + K|2?i(k + K) - if 2^(K)) 



In the case of gravity, using the same approximation as 
in Eq. (|52|l we conclude that 



^w^(k) = -nok^w(k) = AnGpo. 



(B4) 



n=l 
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